function correctionMask=inhomogenityCorrection(image,sigma)
%% Gaussian smoothing of 3D image in fourier domain.
% Correction smoothing radius can be modulated with the sigma variable.
%
% Original script from Prof. Norbert Gretz's Lab in Heidelberg, Germany 
% (website http://www.umm.uni-heidelberg.de/inst/zmf/index_e.html).
% Authors of original script: Ivo Wolf and Lothar Schad.

% Modifications by Michael Eager
% Monash Biomedical Imaging, Monash University

%     Copyright © 2012-2013 Ivo Wolf, Lothar Schad, Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.


sz_image = size(image);
if nargin < 2
    sigma = 3./8;
end
zMaskSize = 5;

% Create 3D filter in freq domain
  fG=fspecial('gaussian', sz_image(1:2), 1/sigma);
  fG=bsxfun(@times,fG,ones([1 1 sz_image(3)]));
  clear hz;
  hz(1,1,:)=fspecial('gaussian', [sz_image(3) 1], 1/sigma);
  hz=bsxfun(@times,hz, ones([sz_image(1) sz_image(2) 1]));
  fG=hz.*fG;
  clear hz;
  fG=ifftshift(fG);
  
  
% do smoothing in fourier domain
  fSm=fftn(image);
  fSm=fSm.*fG;
  clear fG;

  
% back to image domain
  fsmoothed=real(ifftn(fSm));
  clear fSm;

  % normalization factor of fG is not included: not completely clear for
  % even image dimensions => it is easier to do it afterwards.
  correctionMask=fsmoothed/mean(fsmoothed(:));
  clear fsmoothed;

end

